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Abstract. We consider the nonequilibrium dynamics of an interacting spin-| 
fermion gas in a one-dimensional optical lattice after switching off the confining 
potential. In particular, we study the creation and the time evolution of spatially 
separated, spin-entangled fermionic pairs. The time-dependent density-matrix 
renormalization group is used to simulate the time evolution and evaluate the 
two-site spin correlation functions, from which the concurrence is calculated. We 
find that the typical distance between entangled fermions depends crucially on 
the onsite interaction strength, and that a time-dependent modulation of the 
tunneling amplitude can enhance the production of spin-entanglement. Moreover, 
we discuss the prospects of experimentally observing these phenomena using spin- 
dependent single-site detection. 
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1. Introduction 

1.1. Motivation 

The tremendous experimental progress with ultra cold atoms in optical lattices makes 
them a unique playground for studying nonequilibrium dynamics of many-body 
systems. The remarkable features of these systems are the precise and dynamical 
control of the interparticle interaction and external trapping potentials, as well as 
long coherence times [lj. Different nonequilibrium initial states can be generated by 
a quantum quench (for a recent review see |2J): a parameter in the Hamiltonian 
is suddenly changed such that the system, initially in the ground state of the 
Hamiltonian, is afterwards in an exited state of the new Hamiltonian. Quenches have 
been experimentally realized, for instance, in the interaction strength [31 HI [5] and in 
the additional trapping potential by either switching it off [6] or displacing its center 
[8j [9] . This led to observations such as the collapse and revival of the coherence in 
a Bose-Einstein condensate after a quench from the superfluid to the Mott insulating 
regime [3j. 

New detection schemes [lOj [TTJ [12] provide the possibility of observing the many- 
body state at the single-site and single-atom level and make these systems even 
more suitable for studying the dynamics of (spatial) correlations in nonequilibrium 
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situations. These methods have already been used to monitor the propagation of quasi- 
particle pairs [T3] and a single spin impurity [15J in a bosonic gas. They have also 
inspired theoretical work on many-body dynamics subject to observations, treating 
issues such as entanglement growth in a bosonic system |16| . destructive single-site 
measurements [T7] , and Quantum Zeno effect physics for spin [18] and expansion [T9] 
dynamics. 

In the present article, we study a different aspect of these systems: the creation 
and dynamics of spin-entanglement during the expansion of a strongly interacting, 
spin-balanced fermionic gas in an optical lattice. We find that the expansion out of an 
initial cluster of fermions can automatically generate long-range spin-entanglement. 

Expansion dynamics of interacting fermions in an optical lattice has been realized 
recently in a three-dimensional optical lattice |6j. In that experiment, the authors 
observed a bimodal expansion with a ballistic and diffusive part and were able to 
explain its anomalous behaviour using a Boltzmann-based approach. In addition, 
expansion dynamics of this kind has been studied numerically in one dimension, going 
beyond a kinetic description. These studies revealed the dependence of the expansion 
velocity [20j I2T] . the momentum distribution function, and the spin and density 
structure factors [55] on the onsite interactions and the initial filling. Furthermore, the 
effects of the different quench scenarios [23] and of the gravitational field [53] on the 
time-evolution of the density distribution have been discussed. The sudden expansion 
of a spin-imbalanced fermionic gas has been recently considered for observing Fulde- 
Ferrell-Larkin-Ovchinnikov correlations [25 , 26| [27]. 

In general, the dynamics is crucially affected by the difference in behaviour 
between a single fermion and that of a doublon (i.e., a pair of fermions at the same 
lattice site). For large interaction strengths doublons are very stable against decay 
into fermions (analogously to the repulsively bound boson pairs |28j) and move slowly. 
These properties can lead to the condensation of doublons [29] and a decrease of the 
entropy [30] in the center of the cloud during the expansion. In the present work, we 
will show how spin-entanglement is generated by the decay of a doublon into single 
fermions. 

The article is organized as follows: First, we will describe in detail the expansion 
protocol, starting from a finitely extended band insulating state, i.e., a cluster of 
spin-singlet doublons. Moreover, we give details about the numerical time-dependent 
density matrix renormalization group (tDMRG) simulation. Before studying spin- 
physics, we first discuss correlations in the density of the expanding cloud (section |2|, 
where the influence of interactions is already significant. We find that large onsite 
interactions lead to positive correlations between certain lattice sites. The dynamics 
of spin-entangled pairs, the main focus of our work, is then presented in section [3j 
First of all, we relate the concurrence to the spin-spin correlation functions. Then we 
discuss the propagation of spin-entangled pairs in the cloud. Furthermore, we compare 
the cumulative "amount" of spin-entanglement in different regions of the cloud and for 
various interaction strengths. We find that spin-entanglement between remote lattice 
sites is only found for large interaction strength, while for small onsite-interactions 
there is entanglement preferentially only between nearby sites. Moreover, the Pauli- 
blocked core of the cluster favours both partners of a spin-entangled pair to be emitted 
into the same direction when compared to the decay of a single doublon, where 
they almost always are emitted into different directions. In addition, we discuss the 
expansion under the action of a time-dependent (modulated) tunneling amplitude. We 
find that the modulation can enhance the production of spin-entanglement. Finally, we 
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Figure 1. Schematic of the expansion from a band insulating initial state into 
an empty one-dimensional lattice, (a) Due to the discreteness of the lattice the 
velocity of the fermions is bounded. Thus, the fermions emitted from a doublon 
move within a light cone (shown for the leftmost and rightmost doublon). (b,c) 
Examples of fermion configurations at different evolution times, where only the 
rightmost doublon decayed into single fermions and the other doublons remained 
at the original position, (b) A spin up (down) excitation can propagate through 
the band insulating cluster as a spin down (up) hole (time evolution from bottom 
to top), (c) In this case, both fermions have been emitted to the same direction. 



discuss a possible experimental setup for observing this spin-entanglement dynamics 
in the near future using spin-dependent single-site detection (section [3]). 

1.2. Model 

In this paper, we consider an ultracold gas of fermionic atoms loaded into a one- 
dimensional optical lattice. The atoms can be prepared in two different hyperfine 
states, which we label t an d i, and fermions of different "spin" interact via s-wave 
scattering. This fermionic mixture represents a realization of the standard fermionic 
Hubbard Hamiltonian [3U 132] 

L-l L 

i=l a=t,4- i=1 

The first term describes the tunneling of fermions between adjacent lattices sites with 
amplitude J. The second term encodes the effective onsite interaction energy U of 
fermions with different hyperfine states, which can be controlled using a Feshbach 
resonance. The particle number operator is n^ a = cj a c^ a , with the creation 

(annihilation) operator cj a (c^ a ) satisfying the fermionic commutation relations. 

In the following, we focus on the expansion from a band-insulating state, i.e., the 
sites in the center of the lattice are doubly occupied. Experimentally, this is achieved 
by using an additional trapping potential to confine the fermions to a central region 
of the optical lattice [B]. At time t = the trapping potential is switched off and the 
fermions expand into the empty lattice sites, as depicted in figure [TJa). 

The time evolution of the average fermion density and the average doublon density 
has already been studied for this expansion protocol, by numerical simulations for 
ID lattices [20]. For large onsite interactions strengths U/J > 4 two wave fronts 
with different velocities are visible in the time-dependent density profile, while there 
is a single wave front for small onsite interaction strengths. It turns out that the 
expansion can be basically understood as a mixture of propagating single fermions 
and doublons, see also [21] . For small onsite interaction strengths the initial state 
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quickly decays into single fermions, which move ballistically through the lattice. Due 
to the cosine dispersion relation of a single particle in the lattice, e& = — 2 J cos (k) with 
wavenumber k G (— 7r,7r] , its maximal velocity is given by |v ma:E | = 2 J. This fact 
leads to light cones in the density distribution, as indicated in figure [TJa). For large 
onsite interaction U/J > 4, only a small fraction of doublons decay into fermions, 
which move away ballistically. As the effective tunneling amplitude of a doublon is 
2J 2 //7, they initially remain in the central region and then move through the lattice 
on a much larger time scale. 

However, the propagation of the single fermions is highly correlated, as they are 
always created as spin-singlet pair when a doublon decays. As sketched in figures [T](b) 
and (c), the two particles of such a spin-entangled pair may be either detected on the 
same side of the initial cluster or on different sides. 

We briefly point out that the dynamics of the spin-entanglement as well as that 
of the density-density correlation is invariant under the change of the sign in the 
interaction strength U. This is a direct consequence of a transformation property of 
the (spin- or density-) correlators and of the initial state under time reversal and tt- 
boost (translation of all momenta by it). This dynamical U i->- — U symmetry in the 
Hubbard model is discussed in more detail in [53]. 

1.3. Numerical simulation 

For the numerical evaluation of the nonequilibrium time evolution, we employ the time- 
dependent density-matrix renormalization group (tDMRG) method [34 , 35 , 36 , 37, 38J. 
The initial state is a cluster consisting of 10 doublons, which are located in the center 
of a lattice of size L — 100 with open boundary conditions. Our tDMRG simulation 
uses a Krylov subspace method with time step JSt = 0.125 and the discarded weight 
is set to either 10 -5 or 10 -6 , depending on the interaction strength. For all evolution 
times shown in the figures, the density remains negligible at the boundaries of the 
lattice. We also verified that the features presented here does not dependent on the 
precise position of the cluster in the lattice nor change when the truncation error is 
modified within the given range. 

For non-interacting particles, we have used the exact expressions for the time- 
dependent correlation functions and concurrence, as given in appendix A. 

2. Density-density correlation 

Before turning to the spin-entanglement (in the subsequent section), we will first 
study the density-density correlation function Dij(t) = (hi(t)hj{t)) — (hi(t))(hj(t)} , 
with hi = hi^ + fii^. For the expansion, sketched in figure [Tj the density at different 
lattice sites is expected to be correlated for several reasons, in particular since the 
fermions are always created in pairs out of doublons. 

Note that the single-site detection of density correlations (more precisely the 
parity correlations) has been very recently used to study the quasi-particle propagation 
in a commensurate ultracold bosonic gas after an interaction quench [T3l [T4] . It has 
also been employed to study the role of onsite interaction in a bosonic two-body 
quantum walk, experimentally realized in a nonlinear optical waveguide lattice [39J. 

As a point of reference, we first consider the density-density correlations for the 
initial state consisting of a single doublon. The fluctuation in the density at a lattice 
site z, Da(t), is the variance of the occupation number hi. For the expansion from a 
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Figure 2. Off-diagonal density correlations Dij^i(t) = (hi(t)hj^i(t)) — 
(hi(t)) {fij^iit)) for the expansion from an initial state consisting of a single 
doublon (a,b) and a cluster of ten doublons (c-h). For better visibility we display 

1 /3 

the values of as a colour scale and set Da to zero. (c,f) Initially all fermions 

but the outermost ones are Pauli-blocked. This leads to non-vanishing density- 
density correlations only close to the edge of the initial configuration for small 
evolution times. (a,c-e) For vanishing onsite interaction U, the density-density 
correlation equals the two-site spin-z correlation and is always nonpositive, cf. 
equations ( |A.5| ) and ( |A.7| ). (b,f-h) In the interacting case, the density-density 
correlation can be positive, (b) For a single doublon, Dij^i(t) > only for i 
and j at different sides of the initially occupied lattice site. Thus, the doublon 
decays primarily into fermions moving in opposite directions, (f-h) For a cluster 
of doublons, positive density correlations are also found for lattice sites i and j 
that are both located to the left or to the right of the initial cluster position, or 
within the initial cluster (indicated by the dashed square). 



doublon, it is maximal for those lattice sites located at the edge of the single fermion 
light cone (doublon light cone) in the noninteracting (strongly interacting) case. In the 
following we focus on the off-diagonal density-density correlations D^^(t), which are 
shown in figures |2|a) and (b). Most importantly, we find that the onsite interaction 
U leads to a positive correlation of those fermions propagating with almost maximal 
velocity \2J\ in opposite directions, see figure |2jb). On the other hand, the density- 
density correlation assumes large negative values between those lattice sites in the 
center that are most likely occupied by the doublon (and not by single fermions) 
[central square region in figure [2jb)]. In contrast, the off-diagonal density-density 
correlations are always nonpositive for clusters of noninteracting fermions, as detailed 
in appendix A, see also figures |2ja) and (c)-(e). 

The bunching effect in the interacting case can be understood by writing the 
motion of the two fermions in relative and center of mass coordinates, r and R, 
respectively (see appendix B for details of the discussion). Note that we assume 
an infinite lattice for this argument, which is compatible with the simulation as the 
boundary conditions do not play a role for the evolution times considered here. The 
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center of mass motion is a plane wave with total wavenumber K = (k\ + fe) mod27r, 
where k\ and ki are the asymptotic wave numbers of the single fermions. The relative 
motion is described by a if -dependent Hamiltonian. For [/ / 0, the eigenstates of 
this Hamiltonian are one bound state and scattering states. The probability for the 
doublon to decay into one specific scattering state with wavefunction ipK,k(r) [where 
k = (ki — &2)/2 is the relative wavenumber] is given by the modulus squared of that 
wavefunct ion's probability amplitude at r = 0: |^k,/c(0)| 2 . This decay probability is, 
up to an overall normalization constant [see equation ( |B.5| )]: 

fc(0)| 2 oc [1 + U 2 / (16J 2 cos 2 (K/2) sin 2 (/c))] . (2) 

For small onsite interactions \U/4J\ <C 1, the decay probability is almost the same for 
the different scattering states, except for K « tt or k ~ where it drops to zero. In the 
strongly interacting case, it exhibits a pronounced maximum for K = and k = ±7r/2, 
that is k\ = ±7r/2 and ki — T^/Z. In other words, for large onsite interaction the 
doublon decays primarily into two fermions moving in opposite directions with velocity 
close to |2J|. 

Next we study the density correlations for the expansion from a cluster of several 
doublons into an empty lattice. The results are displayed in figures [2jc)-(h). Just 
as for the single doublon, the off-diagonal density correlation has regions of positive 
values in the presence of onsite interaction. However, in contrast to the case of a 
single doublon, the existence of other doublons now leads to positive correlations also 
at lattice sites located on the same side of the initial cluster, see figures [2jf)-(h). 
Positive correlations between sites on different sides of the cluster are only observed 
for evolution times larger than the time a hole takes to propagate through the cluster, 
cf. figure [2^h). 

The results suggest that the initially Pauli blocked core of the cluster leads to an 
enhanced decay of the outermost doublons into single fermions moving away in the 
same direction. However, the alternative case, where the edge doublon decays into a 
fermion and a hole moving into the opposite direction close to the maximal velocity 
\2J\ also leads to positive correlations. Further simulations show that positive density 
correlations between lattice sites on the same side of the initial cluster positions are 
found for clusters of about four and more doublons and U/J > 2. Note that as the 
single fermions are created by the decay of an edge doublon, the situation is different 
from a free single fermion or fermionic wave packet approaching a cluster of doublons. 
In the latter situations, the fermion is almost perfectly transmitted through the cluster 
(band insulating state) in the limit of large onsite interaction [4Q| I4T) . 

In the next section, we discuss the spin-entanglement. This will reveal that the 
positive density correlations indeed stem from singlet pairs, which become delocalized 
by the decay of a doublon. 

3. Spin-entanglement between different lattice sites within the expanding 
cloud 

In this section, we discuss the entanglement between fermions located at different 
lattice sites by means of the concurrence [32 j. First, the relation between the 
concurrence and the spin-spin correlations is established. Then we use tDMRG 
simulations to find the regions where fermions are likely and unlikely to be entangled 
during the expansion. We examine the role of the onsite interaction and the core of 
the cluster on the entanglement dynamics. 
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3.1. Reduced density matrix and concurrence of two fermions 

In this paragraph we derive the reduced density matrix and the concurrence for 
fermions located at different lattice sites within the expanding cloud (see [43] for 
a related discussion in the context of solid state physics). Given the full many- 
body wavefunction the reduced density matrix for the lattice sites i and j 
is Pij(t) = Tr S ites^i,j{|^ r (^))(^ r (^)|}, where all other lattice sites have been traced out. 
Whenever two fermions are situated at the same lattice site they form a spin-singlet 
pair as their spatial degrees of freedom are symmetric under particle exchange. Thus, 
entanglement in the spin degree of freedom between fermions at two different lattice 
sites can only occur if each lattice site is occupied by a single fermion. Projecting 
Pij(t) onto those states yields 

Here, nf = n^ + — 2n^n^ is the single fermion number operator at site z, which 
projects onto a subspace with exactly one fermion on that site. The normalization 
factor in the denominator, Tr{p^(£)nf n|} = (nf (t)n^(t)), is the probability of finding 
at time t a single fermion at each of the lattice sites i and j. The state described 
by pfj(t) is the state obtained after a successful projective measurement. The 
reduced density matrix pfj(t) is equivalent to a two-qubit density matrix, which can 

be expressed in the form p = ^ =0 A a/3 <r^ (8) o - ^, where a 1 ' 2 ' 3 are the Pauli 
matrices, a = 1, and the factors A a/3 are determined by the correlation functions 
^a(3 _ |(a^(j^) [44J. Consequently, the reduced density matrix of single fermions 
at lattice sites i and j can be written as 

w-jmrni.®^ (4) 

where Sj ,2jS = \ J2 a ,b=t i ^l,a(^" 1 ' 2 ' 3 ) a «A > & is the x-, y-, and z-component of the 
spin operator and Sf is, for compactness, defined as half the single fermion number 
operator, := ^nf . Note that (5f (t)Sj(i)) is calculated using the full (unprojected) 
wavefunction since states with vacancies or doublons at site i or j do not contribute 
to the expectation value. 

Symmetries of the initial state and the Hamiltonian can simplify the form of 
the reduced density matrix, such that only a few correlation functions are needed to 
determine pfj(t). As detailed now, the reduced density matrix pfj(t) depends only on 
(Sf(t)Sj(t))/(hf(t)hj(t)) for the cluster initial state shown in figure [lj The Hubbard 
Hamiltonian preserves the spin-dependent particle number, i.e., [H,N^^ = 0, 
N-\-,l = ^2f = i nt,l- Given an initial state with fixed number of spin- up and spin-down 
fermions, which is the usually case in cold atom experiments, the time-dependent 
expectation values of operators that do change the spin-dependent particle number 
vanish. This yields (nf (t) Sf v \t)) = 0, (S?' v \t)Slf(t)) = 0, and (Sf(t)S?(t)) - 

(l%(t)SV(t)) ± i[(Sf(t)SV(t)) + (S%(t)S?(t))] = 0. The latter condition comes from 
creating two spin-down (spin-up) fermions while destroying two spin-up (spin-down) 
fermions at lattice sites i and j and can also be written as (Sf(t)Sj(t)) = (Sf(i)Sj(i)) 
and (Sf(t)S](t)) = -(Sf(t)S?(t)). Moreover, the Hamiltonian Q is fully rotationally 
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invariant, [H, Yli=i S^ v ' z ] = 0. If the initial state is rotationally invariant, for instance 
a cluster of doublons, then the many-body state remains SU(2) spin symmetric during 
the time evolution. It follows that (S?(t)S*(t)) = (Sf(t)SJ(t)) = (S% (t) S? (t)) , 

(n|(t)5|(t)) = (ni(t)S?>v(t)) = 0, and (5f(t)5|(t)) = (S^(t)S*(t)) = 0. In 
summary, the reduced density matrix for the expansion from a cluster of doublons 
reads 

Pi * {t) = 4 Wt)W) ® °* +a * +a * ® a i ] 

I f 1 Q &'^h s )(S I (5) 

where |S^) = ^(| tiij) — I >Utj)) 18 the singlet state and |T™) is the triplet state 
with the m denoting the spin projection in z-direction. 

Instead of (S?(t)Sj(t)) one could in principle evaluate the spin-spin correlation 
in any other direction. However, for cold atomic gases in optical lattices the 
correlation function (Sf(t)Sj(i)) = — ^,^(£)][^j,tW — seems to be 

experimentally most realistic to access as it could be obtained from snapshots of the 
spin-dependent single-site detection of the particle number. 

The spin-entanglement between single fermions can be derived from the reduced 
density matrix. In this work, we use the concurrence C{p) [42J to quantify the 
entanglement. Given the time-reversed density matrix p = <j y (g) <j y p*a y ^a^, with 
the complex conjugation p* taken in the standard basis {| tt)? I tl)? I it)? I -14)} > the 
concurrence is defined by C(p) = max{0, V%~ — V^2~ — V^I}? where the A^'s are 
the eigenvalues of pp in descending order. In our case, the concurrence of the reduced 
density matrix ([5| is given by 

r m Jn 1 Jll^lMlX ,a 
Ci jit) = max < 0, 6 — > . (6) 

' A ' \ 2 (nj(t)h°(t)) f V ! 

Spin-spin correlations (Sf(t)Sj(t))/(hf(t)hj(t)} > —1/12 result in a vanishing 

concurrence. The concurrence approaches 1 as (Sf(t)Sj(t))/(hf(t)hj(t)) \ —1/4, 
i.e., when the two fermions detected at lattice sites i and j always have opposite spin. 

The probability that single fermions at lattice sites i and j form a spin-singlet 
pair can be directly read off the reduced density matrix 

ff = Tr [pum^\] = \ - 3 ||ff|t- ( 7 ) 

Each of the spin triplet states is measured with the probability 1/4 + 

(S!{t)SI{t))/(nl{t)n°{t)). 



3.2. Time evolution of the spin- entanglement within the expanding cloud 

Spin-entanglement between different lattice sites, which we denote by i and j, requires 
that both sites are singly occupied as discussed in section |3.1| The probability for 
this is (nf(t)hj(t)), with the single fermion number operator already defined above 
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Figure 3. Density-density correlation of single fermions (nfi^hj,^)) (i.e. 
excluding doublons) for the expansion from an initial cluster of doublons in 
one dimension (the initial position of the ten doublons is indicated by the 
dashed square), (a-c) Without onsite interactions, U/J = 0, the fermions move 
ballistically through the lattice. This is reflected by the fourfold symmetry of the 
correlation matrix shown here, (d-i) When increasing U/J, fermions are created 
only rarely by the decay of a doublon at the edge of the cluster. These two 
fermions move within the same light cone, cf. figure [I] The light cone leads to 
a square shape of the correlation function having its center at the edge of the 
initial cluster, shown by the dotted square in panel (i). The dotted lines indicate 
pairs of coordinates corresponding to fermions emitted with opposite velocities by 
a doublon at the edge of the cloud. Moreover, an increased correlation between 
nearest neighbor lattice sites is found for larger evolution times, see (f) and (i). 



(nf = hi^ + hi^ — 2fii^fii^). Figure [3] shows the numerical results for (nffyfij^fj;)) 
for different evolution times Jt and onsite interaction strengths U/J. 

Shortly after the quench, single fermions are created at the edges of the cluster, 
see figures [3^a),(d), and (g). For the noninteracting system, the fermions escape 
the cluster successively (i.e., starting at the edges, and finally from the center of 
the cluster). They move ballistically and the correlation function displays a fourfold 
symmetric structure, figures [3^ b) and (c). In the interacting case, in contrast, the 
decay of a doublon into fermions is heavily suppressed. This results in a correlation 
function (hf(t)hj(t)) that has its main contributions within two square regions given 
by the light cones of fermions emitted by the outermost doublons of the cluster, see 
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figures [3jf) and (i). Within the light cones, single fermions are more likely to be found 
at lattice sites corresponding to a motion with almost maximal velocity |v ma;E | = 2 J 
into opposite direction. The density correlation between single fermions attains its 
largest values for nearest-neighbor lattice sites within the cloud, cf. figures [3](e),(f),(h), 
and (i). This may be due to virtual transitions between two configurations: either a 
doublon with a neighboring vacancy, or a state of two adjacent single fermions. Note 
that the mean total number of single fermions, X^=i(^f)> approaches within a few Jt 
an almost constant value, cf. area under the curves for Jt = 3.5,7.5 and U/J = 3,6 
in the last column of figure [3] That means only the edges of the cluster evaporate 
for large interactions, releasing a finite number of single fermions. Moreover, (nf (£)) 
becomes relatively flat in the center of the cloud for larger evolution times Jt. 

Let us now consider the spatial distribution of spin-entangled fermions during the 
expansion. For this purpose, the concurrence between two fermions at different lattice 
sites is numerically evaluated using equation (|6|. Figure [4] shows the concurrence for 
any pair of lattice sites i and j, at different times Jt and onsite interaction strengths 
U/J. 

For the expansion of noninteracting fermions, figures [4^a) -(c), the concurrence is 
finite only for nearby lattice sites. It is almost 1 within the outermost wings of the 
expanding cloud. This can be physically understood the following way: Due to the 
Pauli principle, fermions with the same spin become spatially antibunched during the 
expansion. Thus, fermions at neighboring sites are more likely to have opposite spin, 
cf. figures |2jc) -(e). In addition, the outermost region of the cloud lies only within 
the light cones of the doublons close to the edge of the initial cluster. Two fermions 
detected in this region are almost certaitly emitted by the same edge doublon and in 
consequence have a high probability to be spin-entangled. 

When increasing the onsite interaction U, spin-entangled pairs are formed on 
remote lattice sites, too, cf. figures |4^d)-(i). Indeed, spin-entanglement is found 
between lattice sites within and outside the initial cluster position, figures [4^e) and 
(h), as well as on different sides of the initial cluster position, figures [4jf) and (i). 

The figures are a fingerprint of the creation of a counter propagating hole and 
single fermion by the decay of an edge doublon. When the hole moves through the 
cluster, the spin-entanglement with the single fermion outside the cluster is swapped 
sequentially from one fermion to the next in the cluster. In this way fermions become 
entangled that have never been on the same lattice site and have never directly 
interacted with each other. At the end of this process, a spin-singlet pair is created 
with a fermion at each side of the cluster. The concurrence for two fermions on 
different sides of the initial cluster position increases with the interaction strength, 
compare figures gf) andgi). This can be understood by the suppression of the decay 
probability of doublons with increasing interaction strengths. For larger interaction 
strength, a hole is more likely to cross the cluster without being disturbed by another 
hole. 

For nearest-neighbor sites, which have a relatively large probability to be 
simultaneously singly occupied [see figures [3^ f) and (i)], we observe a concurrence 
close to 1. This implies the creation of vacant lattice sites within the cloud during the 
expansion. The singlet pairs on nearest-neighbor sites come from virtual transitions 
between a doublon with neighboring vacancy and a state of two single fermions. We 
verified this behaviour in addition by numerically creating an ensemble of snapshots 
for the distribution of fermions as described in |19| . 

We emphasize that for inhomogeneous systems, such as the one discussed in 
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Figure 4. Concurrence Cij(t) of two (single) fermions located at lattices sites i 
and j after different expansion times Jt, see equation |6j. Initially, ten doublons 
are located at the sites indicated by the dashed square. We emphasize that dark 
red indicates strictly zero concurrence (also see the cuts displayed to the right). 
The black colour code is used whenever the probability of finding fermions at 
sites i and j is too small, {nf(t)hj(t)) < 10 -5 . In those cases, the concurrence is 
not computed as it would become susceptible to numerical inaccuracies. (a,d,g) 
For short evolution times, here Jt = 1, single fermions are mainly created at the 
edge of the cluster since the central fermions are initially Pauli blocked. Fermions 
close to the same edge of the cloud are likely to be entangled as they are likely to 
originate from the same doublon. Fermions at opposite edges are not entangled 
since they can not have been emitted from the same doublon. (b,c) At larger 
evolution times and without onsite interactions, the concurrence is for most 
sites i and j. However, it is close to 1 when two fermions are located at the 
outermost part of the cloud. (e,f,h,i) By increasing the onsite interaction U/J, 
entanglement of fermions across the cluster becomes possible. In this case, the 
concurrence is highest when the fermions result from the decay of a doublon at 
the edge and escape with opposite velocity, indicated by the dotted line [see also 
cut through panel (i)]. Moreover, the concurrence of fermions at neighboring sites 
becomes almost 1. Within the central region, approximately given by the overlap 
of the light cones shown in figure^ the concurrence remains relatively small. 
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Figure 5. Spin-entanglement of a single lattice site i with all the other sites. The 
panels show the time evolution of Ctot,i(t), given by equation |8j, for different 
interaction strengths U /J. (a) In the absence of interaction, it is primarily the sites 
at the edge of the cloud which are spin-entangled. (b,c) For finite interactions, we 
observe the following: At larger times Ct t,i(t) is nearly homogeneous in a central 
region within the doublon light cones [dotted lines in (c)]. Spin-entanglement for 
lattice sites removed from this central region is finite at locations corresponding 
to trajectories of fermions, which have dissolved from the edges of the cluster. For 
the left edge the single fermion light cone is shown as dashed lines. 



this article, the spatially resolved measurement of two-point correlations provides 
more information about the dynamics than structure factors, such as Sk oc 
J2im e~ %k ( l ~ m " > (§i S^j}. While the latter could be used to determine the average spin- 
spin correlation of fermions at fixed distance, it contains no information where these 
spin-correlated pairs are located in the cloud. 

3.3. Summed concurrences 

Above, we examined the spin-entanglement between two lattice sites for fixed time 
points. In this subsection we aim to quantify the spin-entanglement for entire regions 
in the lattice. In particular, we address the questions: Are there sites which share 
more spin-entanglement with the rest than other lattice sites? How does the spin- 
entanglement in different regions built up as function of time? Which locations are 
most entangled in the weakly and strongly interacting case? How does the size of the 
cluster affect these results? 

In the following we discuss the amount of pairwise spin-entanglement of a lattice 
site or a region in the lattice in terms of the summed concurrence. We define it as the 
sum over the concurrences, Cij(t), which are weighted by the probability of detecting 
a single fermion at both lattice sites, (hf(t)hj(t)). The weights are introduced to 
accommodate for the possibility of vacant or doubly occupied lattice sites, in contrast 
to the summed concurrence used in spin systems [45J . For a system consisting of spin- 
singlet pairs whose wavefunctions do not overlap, i.e., Cij(t) is either zero or one for 
all sites i and j, the summed concurrence equals the average number of delocalized 
spin-singlet pairs. Note that the summed concurrence is by no means a measure for the 
total entanglement of the system. It neither includes multipartite entanglement nor 
entanglement in the occupation numbers. For a detailed discussion on entanglement 
in many-body systems we refer to [46J , see also [47j [48l [49j [16] for recent proposals on 
detecting entanglement in cold atom systems. 

Let us consider the spin-entanglement of a site with all the other lattice sites. 
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The summed concurrence for site i is defined by 

CtatM = y E(nmn' s (t))C itj (t). (8) 

The time evolution of Ctot,i(t) is shown in figure [5] for different interaction strengths 
U/J. For noninteracting fermions [figure |5ja)], lattice sites close to the edge of the 
cloud display the strongest entanglement, while sites in the rest of the cloud are hardly 
entangled. For increasing interaction strengths a central region with almost uniform 
Ctot,i(t) builds up during the evolution, see figures [5^b) and (c). For large U/J, we 
find that Ctot,i(t) approaches the expectation value (nf). This turns out to be related 
to the fact that in this case the probability of having two doublons decay is negligible, 
and the contribution comes almost entirely from the decay of a single doublon. 

In figure Qi) it is apparent that onsite interactions can lead to spin entanglement 
across the expanding cluster. In the following we compare the summed concurrences 
of lattice sites on the same side and on different sides of the initial cluster location, 
C ss and Cdsi respectively. They are defined by 

C ss (t)= J2 (nmh^Qjit) (9) 

i+l<j<l V i—l>j>r 

C ds {t)= E (nt(t)hp))C^(t), (10) 

i<l A j>r 

where I and r denote the leftmost and rightmost occupied lattice sites of the initial 
state. Note that nearest-neighbor lattice sites are excluded from C ss (t). That 
means we do not take into account contributions from virtual transitions of doublons 
(decaying virtually into two adjacent fermions) that move away from the cluster initial 
position. In addition, we consider the summed concurrence of sites at fixed distance 
d, Cd(t) = Y^i(™i(t)™i+d(t))Ci,i+d(t), and the summed concurrence of all sites, 
Ctot(t) — ^2i < j(^ L i(t) r ^j(t))Cij(t). The time evolution of these summed concurrences 
is shown in figure [6j 

For noninteracting fermions and small evolution times, the total summed 
concurrence is dominated by contributions from nearest-neighbor sites. For larger 
times more and more spin-entanglement is transferred to fermions found on the same 
side of the initial cluster position (C ss ), see figure|6ja). The spin-entanglement remains 
relevant only for small distances, reflected in Cg(t) = for all simulated times, cf. 
figure [6^d) . 

By contrast, in the interacting case, spin-entanglement is generated via fermions 
propagating away on different sides of the cluster. This is seen as a finite value of Cds (t) 
for times Jt > 5, which is the time a hole needs to propagate through the cluster, cf. 
figures [6^b) and (c). The total summed concurrence and concurrence between nearest- 
neighbor sites quickly settle into a damped oscillation around a constant value. The 
time evolution of the summed concurrence of sites at distance d > 2 is displayed in 
figures |6je) and (f). Cd(t) remains zero for small times, peaks at times Jt slightly 
exceeding d/4, and decreases for larger times. This shows that spin-entangled pairs 
mainly propagate at almost maximal (relative) velocity 4 J. For large U/J, Cd(t) is 
approximately uniform for distances d < AJt and roughly decays as (Jt) -1 for the 
simulated times. Note that C\(t) plays a special role. Its main contribution does not 
stem from "free" fermions. Rather, it is generated by a doublon virtually dissolving 
into adjacent fermions. 
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Figure 6. Time evolution of spin-entanglement in different areas of the expanding 
cloud, (a-c) We show: the summed concurrence Ctot of all lattice sites, the 
summed concurrence C\ between nearest-neighbor sites, as well as the summed 
concurrences of all lattice sites at same side (C ss ) and at different sides (Cd s ) of 
the cloud (see the main text for the definitions). Nonvanishing Cd s is found only 
for finite onsite interaction. Ctot and C\ quickly settle into a damped oscillation 
around a constant value, for U/J>2. The arrows in panel (c) show corresponding 
values of the summed concurrences for the decay of a single doublon, cf. figure |7| 
(d-f) Summed concurrence Cd of lattice sites at distance d > 2. (d) Without 
interactions only close lattice sites are spin-entangled. (e,f) With increasing 
interaction strength U/J, Cd(t) equals zero for times up to Jt « d/4, followed by 
a peak and a decay for larger times. At fixed time Cd is approximately uniform 
for the distances d < AJt. 



Let us finally discuss the impact of the cluster size on the spin-entanglement 
dynamics. For very weak interactions, a larger number of doublons means that more 
delocalized singlet pairs are created shortly after switching off the confining potential. 
For large interaction strengths up to U/J = 40, we simulate the expansion and 
compare the summed concurrences for different cluster sizes, including the case of 
a single doublon. We summarize the results in figure [7J Note that reasonably large 
evolution times ( Jt ~ 20) for the comparison of the summed concurrences are reached 
only for U/J > 6. The summed concurrence of all sites, Ctot: agrees for all considered 
cluster sizes and matches the analytical result for a single doublon, see figure (7^a). 
Apparently, the initially Pauli-blocked core has no effect on the number of created 
single fermions for the considered times. For C ss and C^ s , however, we find a clearly 
different behaviour for a single doublon and a cluster of four and more doublons 
[figure (7^b)]: A cluster is more likely to emit (delocalized) spin-entangled pairs into 
the same direction. 
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Figure 7. Summed concurrences for different sizes of the initial cluster. We 
compute the values at time Jt = 20 (except for U/J = 6, where Jt = 13.5), when 
they are almost constant as function of time, (a) The total summed concurrence, 
Ctot-, shows no dependence on the cluster size for large onsite interaction strengths 
U/J>6. The data agrees with the exact result for a single doublon derived in 
appendix C, Ctot = 1 - [1 + 16J 2 /f/ 2 ]" 1 / 2 , which is shown as solid line, (b) The 
summed concurrences of lattice sites at same side (C ss ) and at different sides (C^ s ) 
of the cloud (see the main text for the definitions) disagrees for a single doublon 
(dots) and a cluster of doublons (squares). Note that clusters of sizes 4,7, and 
10 give similar values. For strong interactions, a cluster prefers the emission of 
(delocalized) singlet pairs into the same direction compared to a single doublon. 
The summed concurrence between nearest-neighbor sites (Ci) approaches Ctot 1 1 
in both cases. Solid lines show analytical results for C\ as well as the contribution 
of scattering states to C ss and Cd s for a single doublon (see appendix C). 



3.4- Expansion with modulated tunneling amplitude 

In the previous section we have seen for the interacting case that total summed 
concurrence, Ctot, approaches a fixed value shortly after switching off the potential, 
via the escape of a few fermions from the cluster edges. Here, we discuss a way of 
"continuously " generating single fermions and enhancing C to t compared to the free 
time evolution. We consider an expansion during which the tunneling amplitude 
is repeatedly varied in time, while the interaction strength is constant. Such 
modulation may be experimentally realized by either varying the laser intensity (the 
tunneling amplitude decreases much faster with increased laser intensity than the 
onsite interaction strength, see, e.g., or by shaking the lattice sinusoidally [50J. 
We find for certain values of the tunneling amplitudes and time intervals between the 
quenches an increased amount of spin-entanglement as shown in figure |U For the 
presented case, the repeated quenches lead to the generation of more and more single 
fermions, see figure [8^a). This results in an enhanced spin-entanglement between 
distant lattice sites, while the spin-entanglement between nearest-neighbor sites is 
suppressed [compare figures |8jb) and|8|c) with figures [6^c) and^f)]. 



4. Remarks on observing the spin-entanglement in experiments 

In the main part of this article we have analyzed the dynamics of the spin-entanglement 



for the expansion from a cluster of doublons. As discussed in section 3.1 the 



concurrence between two lattice sites i and j can be determined by the single 
fermion expectation value (nf (t)h s - (t)) and the spin-spin correlation (Sf(t)Sj(t)) = 
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Figure 8. Expansion with time-dependent tunneling amplitude (modulated in a 
stepwise fashion). The tunneling amplitude is repeatedly switched between the 
values J [time intervals marked by gray background in panel (a)] and J' = J/2 
[white background in panel (a)], while the onsite interaction strength is fixed at 
U I J = 6. (a) This dynamics (dashed line) produces a larger total number of single 
fermions (N s (t)) than the free expansion with amplitude J (solid line). (b,c) Time 
evolution of the summed concurrences. In comparison to the free expansion, cf. 
figures |6jc) and (f), the total summed concurrence as well as the concurrences 
C ss and Cds are enhanced, while the summed concurrence at nearest-neighbor 
sites, Ci, is decreased. This is also seen in the enlarged summed concurrence of 
sites at distances d > 2. 



i([n a (t)-n a (t)][n i>t (t)-n M (t)]). 

Experimentally, both expectation values could be obtained by averaging over 
many snapshots of the spin-dependent single-site fermionic particle number analogous 
to the already implemented single-site detection of bosonic particles [TTJ [12]. Since 
doubly occupied lattice sites do not contribute to these expectation values, it would 
suffice to be able to detect singly occupied sites. Thus, the loss of atom pairs due to 
inelastic light-induced collision during the imaging process, showing up in the bosonic 
case, would not be an issue. 

Measuring correlations of the spin-z component will be sufficient for obtaining 
the concurrence under the following assumptions: (i) the initial state is of the type we 
have described (total singlet, i.e. total spin 0); (ii) the dynamics proceeds according 
to the model Hamiltonian, i.e., a SU(2) symmetric Hamiltonian with no decoherence 
or entanglement with external degrees of freedom. A scenario where the initial state 
violates condition (i) would be a cluster with impurity sites that contain only single 
fermions. If there is only one impurity containing a single fermion, then the problem 
can still be diagnosed because the final snapshot will reveal unequal numbers of spin- 
up and spin-down fermions. In general, however, if these two conditions are in doubt, 
the experiment would have to measure the full spin-spin density matrix of two lattice 
sites, by repeated runs and measurements of spin projections in different directions. 
This could be done by implementing a rotation in spin space before measurements, 
in analogy to the state tomography realized with trapped ions [5T] , While being 
more challenging than measuring the spin-dependent fermion number, such a kind 
of coherent spin control seems to be possible in future experiments. Spin flips at 
single lattice sites have already been shown experimentally for ultracold bosons [52J. 
This setup could be in principle extended to coherent spin control, replacing the Rabi 
frequency sweep by driving a Rabi oscillation [53] . 
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5. Summary and Outlook 

In this article, we have analyzed the creation and time-evolution of spin-entanglement 
during the sudden expansion from a cluster of doublons into an empty lattice. 
Interestingly remote spin-entangled pairs are created for large onsite interaction. In 
addition, an extended cluster favours the emission of the two fermions of a pair into 
the same direction when compared against the decay of a single doublon. Finally, we 
found that a time-dependent modulation of the tunneling amplitude can be used to 
increase the "production" of spin-entangled pairs. Our results provide a starting point 
for studying the spin-entanglement dynamics for more complex initial states or, e.g., 
for spin-imbalanced fermionic gases. 

In the scenario considered here, the spin-entanglement can be extracted from the 
two-site spin-z correlation functions. Thus, it will become experimentally accessible 
once spin-dependent single-site detection has been implemented. 
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Appendix A. Correlation functions and concurrence of noninteracting 
fermions 

In this appendix, we derive analytical formulae for the correlation functions of 
expanding, noninteracting fermions. For vanishing onsite interaction, fermions of 
different species are uncorrelated and each fermion propagates with the free dispersion. 
The time evolution of the annihilation operators is given by 

A?> W = G 3™(t)Cm,a, (A.l) 
m 

with spin index a =t, 4- and the free fermion propagator 

/ n dk 
— exp{i[(m-j)k-2Jcos(k)t]} = V~ m ^_ m (2Jt).(A.2) 

Here, J denotes the Bessel function of the first kind and we assumed an infinite lattice. 
For the previously discussed initial state of localized doublons, the Green's function 
reads 

Sij(t) = (4,a(t)Cj,a(t)) = * J ' _< Jj-m(2Jt)Ji- m (2Jt), (A.3) 

rneO 

where the summation is taken over all initially occupied lattice sites O. For both spin 
species, the density is given by 

W) = G«(t)= J2 J \ 2 i-m\(2Jt)- (A-4) 

rneO 

When the fermions are initially localized at the lattice sites, equal time, normal 
ordered products of operators can be expressed as a Slater determinant of the equal 
time one-particle Green's functions, Gij(t), [54J. The correlation function {Sf{t)Sj{t)) 
can be evaluated by writing it in terms of density-density correlations and making use 
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of (fii^nj^i^t)) 
yields 



Ni{t)Nj{t) - 5 ab |0y(t)| 2 and (n i>a (t)n i>b (i)> = 



2Sab) 



(smsm) 



i 



-\\Q*At)\\ 
NS) [i-Ni{t)}. 



This 

(A.5) 
(A.6) 



Note that the probability of finding two fermions with the same spin at lattice sites i 
and j is by 2\Qij(t)\ 2 smaller than the probability for fermions with antiparallel spin. 
This difference leads to a nonvanishing spin-spin correlation even in the absence of 
onsite interactions. For noninteracting fermions, the spin-spin correlation is related 
to the density-density correlation Dij(t) = (hi(t)hj(t)) — (fii(t))(hj(t)) via 



D ij (t)=A(smsm- 



(A.7) 



Thus, the density-density correlation Dij(t) is smaller or equal to zero for different 
lattice sites i and j, cf. figures |2| a) and (c)-(e). Analogously, the density-density 
correlation of single fermions is calculated using the relations for {ni^ a {t)hj^{t)) given 
above. We find 



-\Gij(t)\ 2 ) 



(A.8) 
(A.9) 



2|0y(*)| 2 



(ht(t)nm) =2M(*)[1-M(*)], 
(nf(t)n^(t)) = 4 [Ni(t) - (WWiW 

'_Nj(t) - (mWj(t) - 

The c oncu rrence CjA t) [defined by equation (|6|] can be evaluated using the expres- 
sions (A.5) and (A.9). It approaches one in the limit (^f , -(t))/(n|(t)n^.(t)) \ 
— j. This is the case for |^-^(f)| 2 Mi{t)Mj^%{t) (note that |^^(t)| 2 can only 
assume values between and 1/2). 



Appendix B. Decay of a doublon into scattering states 

This appendix presents the derivation of the decay probability of a fermionic doublon 
into different scattering states \ipK,k}- In doing so, the wavefunctions of the scattering 
states are calculated mainly adopting the procedure for two-particle states in the Bose 
Hubbard model presented in [55] . 

For one spin-up and one spin-down fermion, and an infinite lattice the 
Hamiltonian can be expressed in the form 

u D = - jJ2i\ ti><ti+i i + 1 ti+i)<ti i + 1 u)(u+i i + 1 u+i)(u i} 

i 

+ uj2\^u)(-ti,u\. (b.i) 

i 

Analogously, we write the two-fermion states in terms of the basis {| ti>ij)}> 
= ■ ^(ti, tii-lj)- For the decay of a doublon, the fermions are always in the 
spin-singlet subspace. Thus, we can write the two-particle wavefunction as ^(ti, -Ij) — 
(Psitil) ' with the spin-singlet wavefunction </? s (a, b) = ^ a ,t^,l — ^a,|^6,t an d a 

symmetric spatial wavefunction, ^(i, j) = Plugging into the Schrodinger 

equation for Hamiltonian ( |B.1| ) yields 

(E - USij) = -J m 1>(i + + 1>(i, 3 ~ 1) 

+ip(i,j + l)\. (B.2) 
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This relation is simplified by introducing center of mass (c.o.m.) and relative 
coordinates, R = (i + j)/2 and r = i — j, respectively. The wavefunction factorizes 
into a plane wave motion of the c.o.m. with total wavenumber K G [— 7r, 7r) and a 
if-dependent relative motion, i.e., = e* XjR ^x(r). The relative motion satisfies 

the recurrence relation 

- Jk [*l> K (r - 1) + ^k(t + 1)] = (E K - US r0 ) il> K (r), (B.3) 
with if-dependent tunneling amplitude Jk = 2J cos(K/2) and energy Ek- For 



vanishing interaction strength U equation (B.3) is solved by plane waves ^K,k(r) = 
e ±ikr w ^h corresponding eigenenergies Ex,k — — 2Jk cos(k). In the interacting case, 
we make a scattering ansatz by writing the wavefunction as a superposition of incoming 
and reflected plane waves, ^K,k(r > 0) = e %kr + ce~ zkr . Here, k is the relative 
wavenumber and ^K,k(r < 0) is determined by the symmetry condition ip(r) = ip(—r). 



The boundary condition at r = in equation ( |B.3[ ) fixes the coefficient c and we obtain 



cos(/cr) + 97 sin(fe|r|) 



(B.4) 



Finally, we compare the decay probability of a doublon for different scattering states 
\ipK,k)- I n doing so, we express |^k,/c(0) | 2 m terms of the average density in the relative 
coordinate, n, which is obtained by averaging \^K,k(r)\ 2 over one period 2n/k. Note 
that n does only depend on the systems size and is independent of fc, K, and U. We 
find that the decay probability into the scattering state \ipK,k) equals 

\^kM°)\ 2 = n [1 + U 2 / (16J 2 cos 2 (K/2) sin 2 (/c))] _1 . (B.5) 
Appendix C. Summed concurrences of a single doublon 

This appendix provides the calculation of the summed concurrences for the expansion 
from a single doublon. The results are depicted in figure [7] For this initial state, 
the fermions form a spin-singlet for all times. Consequently, the concurrence Cij(t) 
[defined by equation ([6])] equals one for all sites i and j, and the summed concurrences 
simplify to sums over the expectation values (nf(t)n^(t)), see section |3.3l 

Let us first consider the total summed concurrence for a single doublon C to t(t) = 
^ i< j(nf(t)hj(t)). This is nothing but the probability of finding the fermions at 
different lattice sites. It can be expressed as 1 — with the doublon survival 

probability P]j(t), i.e., the probability of finding the doublon intact at time t. In the 
long time limit, only fermions in a bound state remain localized close to each other. For 
a singlet pair with finite onsite interaction U ^ and c .o.m . wavenumber K, it exists 
one bound state *p b K (r) [localized solution of equation (B.3), details are given below], 



where r is the relative coordinate. Thus, we obtain Pd(oo) = ^ f^dK \i/j k (0)\ 4 
in the limit of an infinite lattice, where i/j k (0) is the overlap between the doublon 
and the bound state. Analogously, we find Ci(oo) = ^ dK |^(0)| 2 |^(1)| 2 , 
where we have used that C\(t) = J2i(nf(t)hf +1 (t)) equals the probability of finding 
the fermions at time t at nearest-neighbor lattice sites. The bound state can be 
calculated using the exponential ansatz ^ b K {r) = l/y/Njc • oc^ for the wavefunction 



in equation (B.3). This gives a K = U/2J K - s\gn(U/2J K ) ■ y/l + (U/2J K ) 2 and 



M K = (2J K /U) • y/1 + (U/2J K ) 2 , with J K = 2cos(if/2). Inserting this solution into 
the doublon survival probability yields P D (oo) = ^Jl 7T dKJ\T^ 2 = [1 + 16J 2 /U 2 )' 1 ^ 2 . 
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In the limit of infinite times, the total summed concurrence and summed concurrence 
of nearest-neighbor lattice sites are, hence, given by 



Ctot(oo) 
Ci(oo) 



1 



[1 + 16J 2 /U 2 
X - J^dK\a K \ 2 /V } 



1/2 



= 8J 2 /U 2 
AJ 2 /U 2 - 



- 0([J/U} 4 ) 
■0([J/U} 4 ). 



(C.l) 
(C.2) 



In the long-time limit, the summed concurrences of lattice sites at same side 
(C ss (t)) and at different sides (Cd 3 (t)) of initial doublon position can be related to 



the scattering states calculated in appendix B. We denote by C 



(scat) 



and C, 



(scat) 



the contributions to the summed concurrences stemming from the scattering states. 
C { d s s cat) is the sum of the occupation probabilities of scattering states corresponding 
to fermions moving in opposite direction [fci G (0,7r), k<i G (— 7r,0) or k\ G (— 7r, 0), 
^2 G (0, 7r)], and CsT at ^ is the sum of the occupation probabilities of those states with 
fermions moving into the same direction [&i,/£2 G (0,7r) or fci,fc2 G (— 7r, 0)]. Here, 
k\ and &2 denote the asymptotic wavenumbers of the two fermions. The occupation 
probabilities are the absolute square of the overlap between the doublon an d th e 
scattering state, \i/jK,k(0) | 2 - The explicit form of \ipK,k(Q) | 2 is given by equation (B.5). 
Taking the continuum limit we find 



c. 



(scat) 

ds 



(2tt)2 



/ dfci / dfc 2 + / dfci / dk 2 

JO J -TT J -TV JO 



(C.3) 



[1 + U 2 / (16J 2 cos 2 ([fci + k 2 }/2) Bin 2 ([A!i - fc 2 ]/2))] 1 



4^[l/2 + 4/7r 2 ] + 0([J/[/] 4 ), 



(C.4) 



^j(scat) 



wAL dki L dk2+ L dh L 



(2tt)2 

x [l + C/ 2 /(l6J 2 cos 2 ([fci 
4^[l/2-4/7r 2 ] + 0([J/C/] 4 ). 



dk 7 



(C.5) 



fc 2 ]/2) sin 2 ([Ai-fe]/2))]" 



(C.6) 



The summed concurrences C ss (t) and Cd 3 (t) used for the numerical data [equations ([9| 
and ( [To] )] contain additional contributions from bound states. In the definition of 
C ss (t) we excluded nearest-neighbor lattice sites in order to remove most of these 
contributions [note that nearest-neighbor sites do not appear in Cdsif)]. From the 
expression for the bound state given above follows that all other terms are of the 



order 0([J/U] A ). Thus, C ss (t) and C ds (t) agree with C { s s s cat) and C^ Cttc; for large Jt 
and U/J, cf. figure gb). 

In conclusion, we find following relations between the summed concurrences 
for long evolution times and large interaction strengths U/J ^> 1: C\ = 
C ss = [1/4 - 2/7r 2 ]C tot , and C ds = [1/4 + 2/^ 2 ]C £o£ . 
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